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Abstract: We develop a systematic approach to the discovery of parallel iterative 
schemes for solving the shape-from-shading problem on a grid. A standard procedure 
for finding such schemes is outlined, and subsequently used to derive several new ones. 

The shape-from-shading problem is known to be mathematically equivalent to anon- 
linear first-order partial differential equation in surface elevation. To avoid the problems 
inherent in methods used to solve such equations, we follow previous work in reformulating 
the problem as one of finding a surface orientation field that minimizes the integral of the 
brightness error. The calculus of variations is then employed to derive the appropriate 
Euler equations on which iterative schemes can be based. 

The problem of minimizing the integral of the brightness error term is ill posed, since 
it has an infinite number of solutions in terms of surface orientation fields. A previous 
method used a regularization technique to overcome this difficulty, An extra term was 
added to the integral to obtain an approximation to a solution that was as smooth as 
possible. 

We point out here that surface orientation has to obey an integrability constraint if it is 
to correspond to an underlying smooth surface. Regularization methods do not guarantee 
that the surface orientation recovered satisfies this constraint. Consequently, we attempt 
to develop a method that enforces integrability, but fail to find a convergent iterative 
scheme based on the resulting Euler equations. We show, however, that such a scheme 
can be derived if, instead of strictly enforcing the constraint, a penalty term derived 
from the constraint is adopted. This new scheme, while it can be expressed simply and 
elegantly using the surface gradient, unfortunately cannot deal with constraints imposed 
by occluding boundaries. These constraints are crucial if ambiguities in the solution of 
the shape-from-shading problem are to be avoided. 

Different schemes result if one uses different parameters to describe surface orienta¬ 
tion. We derive two new schemes, using unit surface normals, that facilitate the incorpo¬ 
ration of the occluding boundary information. These schemes, while more complex, have 
several advantages over previous ones. 
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1. Introduction 

We begin by reviewing the shape-from-shading problem, its formulation as a minimization 
problem, and the use of the calculus of variations in deriving the partial differential 
equations governing the solution of the minimization problem. 


1.1. Preview 

The first study of the shape-from-shading problem was undertaken by Horn (1970 & 
1975). There, the partial differential equation in surface elevation fundamental to the 
problem was converted to an equivalent set of five ordinary differential equations called 
the characteristic strip equations. Algorithms based directly on numerical solution of the 
discrete approximations of these equations arc inherently sequential in nature and have 
difficulty with unavoidable noise in the image data. 

Later, a method lending itself to parallel solution on a grid was developed by Strat 
(1979) using minimization in the discrete domain. Strat used the gradient to express 
surface orientation and so was unable to deal with occluding boundaries, which are known 
to provide crucial constraint needed to avoid ambiguity in the solution, as shown by Bruss 
(1983). For this reason, another approach, based on the stereographic projection of the 
Gaussian sphere, was explored by Ikeuchi and Horn (1981). The calculus of variations 
was used there for the first time in the analysis of the shape-from-shading problem. Their 
method depended on the use of a regularization term in the functional to be minimized. 

In this paper, we carefully examine the role of the variational calculus in the derivation 
of iterative schemes for shape from shading. Previous methods are discussed in detail, 
and rationalized in terms of the new point of view, where appropriate. The application 
of regularization techniques to well-posed problems is called into question. 

We note in particular that the surface gradient should satisfy an integrability con¬ 
straint. Guided by this observation, we attempt to impose integrability in a strict sense. 
We are, however, unable to derive a convergent iterative scheme based on the appropriate 
Euler equation. We learn that such a scheme may be found if we instead incorporate a 
penalty term based on the integrability constraint. This we demonstrate first using the 
gradient to specify surface orientation, as has been customary. The resulting iterative 
scheme is shown to be related to that developed by Strat. 

As already stated, use of the surface gradient precludes incorporation of the occluding 
boundary information. We overcome this difficulty by taking the novel approach of 
adopting surface-normal vectors directly. This leads to iterative schemes that are more 
complex, but manageable. We finally develop two such schemes that: 

• ensure the result is (at least approximately) integrable, 

• avoid the smoothing introduced by a regularizing term, and 
® permit use of the known normals on the occluding boundary. 

None of the previous methods combined all of these features. 
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1.2. The shape-from-shading problem 


Monochrome images of smoothly carved surfaces with homogeneous reflecting properties 
commonly exhibit a variation in image irradiance, or shading. This is due to the inter¬ 
action of four principal factors: the illumination, the shape of the surface, the reflecting 
characteristics of the material, and the image projection. The shape-from-shading prob¬ 
lem may be regarded as that of extracting the shape information encoded in the irradiance 
data. It therefore entails inversion of the image-forming process. 

Because a number of factors are confounded in irradiance values, the shape depicted in 
an image cannot be determined unless additional information is provided. Of considerable 
utility in this regard has been the reflectance map (Horn, 1977), which specifies the 
radiance of a surface patch as a function of its orientation. The reflectance map can 
be computed from the bidirectional reflectance-distribution function and the light-source 
arrangement (Horn & Sjoberg, 1979). Usually it is more practical to determine the 
reflectance map experimentally, by means of a calibration object of known shape, for 
example. In any case, the reflectance map encodes, inextricably, information about the 
reflecting properties of the surface and the distribution and intensity of the light sources. 

In adopting the reflectance map, we implicitly make the assumption that, for the given 
scene conditions, the radiance emanating from a small surface patch is dependent only 
on the orientation of the patch, and not its position in space. This requires that the light 
sources and the viewer be distant. We also assume that the image is formed by ortho¬ 
graphic image projection, and that the surface has homogeneous reflecting properties 1 . 

Formally, given an image, E, and a reflectance map, R, the shape-from-shading prob¬ 
lem may be regarded as that of recovering a smooth surface, z, satisfying the image 
irradiance equation 

E(z,y) = R(z x (x,y),z y (x,y)) 

over some domain Q of the image. Any given conditions on z on the boundary <9f! of the 
region 0 should also be satisfied. Here z x and z y denote the first partial derivatives of z 
with respect to x and y respectively. Since these derivatives will be used frequently to 
specify surface orientation, it is convenient to introduce the short-hand notation 



and 


q = 


dz 

dy' 


The gradient of the surface z at the point (a :,y) is just (p(x,y),q(x,y)'j. The gradient 
points in the direction of steepest ascent and has magnitude equal to the slope in that 
direction. It is further useful to note that a normal of the once-difFercntiable surface, z, 
at (x,y,z{x,y)y can be written 


n 


-p{x,y), 


-q{x,y),iy 


T 7 

This follows from the fact that (l, 0,p(x, y)) and (0, l,q(x,y)) are tangent vectors and 
that the normal must be parallel to their cross-product. For many purposes one can use 

1 These restrictive assumptions were not exploited in the original work on shape from shading. The 
problem formulation, however, is much easier to comprehend if the reflectance map is introduced, and 
that can be done only if these additional constraints arc imposed. 
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either the surface gradient or the normal to specify surface orientation. Each has its own 
advantages, as we shall see. 

It is customary to choose the direction of projection to be parallel to the 2 -axis. On 
the occluding boundary , the direction of projection is tangent to the surface. That is, the 
normal is orthogonal to a unit vector z, parallel to the 2 -axis. Thus we note that at least 
one of p and q become unbounded on the occluding boundary. 

1.3. Employing the variational calculus 

Suppose we seek, over some domain, a smooth surface satisfying various constraints. It 
is useful to obtain from (the given constraints a non-negative expression that measures 
the departure of a particular surface from a satisfactory solution. We may then search 
for a surface that minimizes the expression. As the value of the expression depends on 
the choice of surface, or function, it is termed a functional. 

The search for a function that minimizes an integral expression is the major concern 
of the calculus of variations (Courant & Hilbert, 1953). Here, we find the valuable result 
that the extrema of functionals must satisfy an associated Euler equation. This equation 
can usually be determined in a straightforward way from the functional. We can, as a 
result, transform our surface-recovery problem from one of minimizing a functional, to one 
of solving one or more partial differential equations. Some of the relevant mathematical 
details are presented in the Appendix of this paper. 

In seeking a surface that best matches the aforementioned constraints, we require a 
global minimum of the corresponding functional. However, Euler equations only specify 
conditions on extremal values. We shall make the strong assumption in this paper that a 
solution to the Euler equation constitutes a global minimum of the functional, satisfying 
the constraints optimally. We shall as a result be deluded if we encounter a surface 
that gives rise to either a local minimum, a local maximum, or an inflexion point in the 
functional, for it too will satisfy the Euler equation 2 3 . The assumption here is difficult to 
avoid, given that we shall be dealing with functionals involving a reflectance map whose 
analytic form may not be known in advance. 

Let us suppose that we obtain from an Euler equation a surface that generates a 
global minimum of the appropriate functional. It may be that the constraints on which 
the functional was originally based are satisfied exactly by this function. However, this 
need not be so. Problems can readily be formulated for which there are no perfect 
solutions. But here we find a very important property of this approach: the surface that 
best matches the constraints will generate a global minimum of the functional. This 
is important to vision problems as they typically involve images that are noisy. Exact 
solutions may not exist in this situation. For example, in the presence of noise, there may 
not be a smooth surface that satisfies the image irradiance equation E(x,y ) = R{p,q) 
exactly. There will, however, be a surface that minimizes the integral of the square of 
the difference between E(x,y) and R(p,q) s . 

2 Because of the use of expressions that arc unbounded above, we shall not encounter solutions generating 
global maxima. 

3 The integral of the square of the difference may have a lower bound that is not attained by any surface. 
In that case, a surface may be found for which the integral is arbitrarily close to that lower bound. 
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It is important to observe that there are typically an infinite number of surfaces 
satisfying the Euler equation. Without further constraint, we do not have a well-posed 
problem. In some cases the original problem includes boundary conditions that, taken 
together with the resulting partial differential equations, lead to a unique solution. In 
the case where the unknown function is unconstrained on the boundary, the calculus of 
variations itself provides so-called natural boundary conditions (see Appendix). 

Care must be taken when formulating the functional to ensure that it provides suf¬ 
ficient constraint, for otherwise there may be an infinite number of solutions even with 
boundary conditions. Such a difficulty may be remedied by the addition of a suitable 
regularization term (Poggio & Torre, 1984). This is discussed in more detail in the Ap¬ 
pendix. 

1.4. A procedure for deriving iterative schemes 

We now consider a way of deriving iterative schemes for recovering surface shape. In the 
event that we seek a surface, z, best satisfying various requirements over 0, we do the 
following: 

(1) Select a functional, F, non-negative over fi, such that 

I{z) = JJ F(x,y,z,..dxdy 

constitutes a measure of the departure of z from an ideal solution. 

(2) Absorb into F any constraint that 0 should satisfy over Cl, using Lagrangian 
multipliers if appropriate. 

(3) If the problem is not well posed as it stands, add a suitable regularization term. 

(4) Find the Euler equation that must be satisfied by the surface z minimizing the 
functional /. 

(5) Determine what boundary conditions are needed to ensure a unique solution. If 
there are no constraints on the function around the boundary <90, determine the 
appropriate natural boundary conditions. 

(6) Develop a discrete approximation of the associated Euler equation, using finite- 
difference methods. 

(7) Design an iterative scheme that converges to the solution of the discrete approxi¬ 
mation of the Euler equation. 

The approach, of course, follows the same pattern if the surface is parameterized in 
a different way 4 . Also, similar results can be obtained by applying the finite-element 
method directly to the functional I. 

As we shall see later, the most difficult step here is typically the discovery of an 
iterative scheme that enables one to recover a solution of the discrete approximation of 

4 The surface may, for example, be parametrized using the gradient (p, q) instead of the surface elevation 
z, for example. 
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the Euler equation. Such a scheme should be efficient, convergent, and preferably lend 
itself to parallel implementation. 

Note that it is better to work with a functional that evaluates to zero for perfect 
solutions. In this way, one is relieved of the onus of showing that there are no unwanted 
surfaces that cause the functional to have a smaller value than that generated by a 
satisfactory solution. An additional advantage of functionals evaluating to zero is that 
one may use them to check how close an iterative scheme is to a solution. This is difficult 
with other functionals, as the minimum value is usually unknown. 


2. Previous work 

Only one shape-from-shading scheme (Ikcuchi & Horn, 1981), prior to this work, has been 
devised by explicit recourse to the calculus of variations. Two other schemes, however, 
(Strat 1979, Smith 1982) can be rationalized by application of the calculus of variations. 
We now examine these three schemes in historical sequence. 

2.1. Strat’s method 

Strat (1979) arrived at his method by application of the standard calculus to the discrete 
domain. We present his analysis here as we wish to show later how it can be related to 
a new scheme we develop using the calculus of variations. Rationalizing Strut’s scheme 
directly in terms of the variational calculus is complicated by the fact that it is based on 
an integral (rather than a differential) integrability term. 

First, let the brightness error at a point ( x,y ) be 

E(x,y) - R(p(x,y),q(x,y)). 

This is the difference between the observed irradiance E(x, y ) and that predicted from the 
estimated gradient (p(x,y), q(x, </)) . In the discrete case, we might consider minimizing 
the total brightness error 5 

n m 

y v ~ r{pij> %)) 

t=l j= 1 

by suitable choice of the gradient at each picture cell in the image 0 . In this vein, then, 
by setting the derivative of the expression with respect to p *./ and q^i equal to zero, we 
obtain, for 1 < k < n and 1 < / < m, the two sets of equations 

ifikl R{Pkli < lkl)') Rp {pkl > Qkl) 0 ) 

(Ekl - R(Pkl,(lkl))Rq{Pkl,<lkl ) = 0 , 

5 For simplicity, we assume a rectangular image region here. There is no loss of generality, however, since 
the sums can be taken over whatever region is desired. 

0 By minimizing an expression containing the sum of the square of the brightness error, we are giving up 
strict enforcement of the image irradiance equation. This seems reasonable, given that neither E nor R 
are known with precision. 
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where R. p and R (j are the partial derivatives of R with respect to p and q respectively. 
These conditions can be trivially satisfied if we choose p lJ and q tJ so that 

Eij• 

Since this equation represents but one constraint on the two unknowns p,y and q tJ , we 
expect that, in general, an infinite number of gradient values will satisfy it, for a particular 
i and j. Many solutions can then be constructed by combining arbitrary choices from 
these sets of possibilities at each picture cell. 

The problem is clearly not well posed as stated. We can, however, make use of the fact 
that the gradients at neighboring points are related. Consider an infinitesimal segment, 
6C, of a curve on the surface. The change in z along the segment is given by 

8z — pSx q 8y, 

where Sx and 8y are the changes in x and y along the segment. The total change in z 
along a curve then is just the integral of (pdx + qdy ). In the case of a closed curve, C, 
this integral should be zero. Thus, if (p(x, y), q(x, y)) is the gradient of a surface z(x,y) 
then 

j> (p( x > y) dx + g(*» y) d v) = o, 

for all closed curves, C , in the region f2 7 . 

Let c denote the spacing between picture cells. Consider an elementary square path, 
with the picture cell (i,j) in the lower left hand corner. If we let i correspond to x and j 
correspond to y, then the integral counter-clockwise around this path can be estimated 
by 

e ij = 2 Pi+ i)i '?»-t-lj “h 9* 1-1 j'+l ~ Pi + lJ+1 ~ Pi,j +1 "" dij+l “ di,j] ■ 

This expression can be obtained by approximating the slope along each of the four sides 
by the average of the slopes at the beginning and end of each side. The result is exactly 
equal to zero when z is quadratic, as can be seen using Taylor series expansion 8 . The 
difference between this expression and the exact loop integral is (perhaps surprisingly) 
of order e 4 . 

On a discrete grid, we wish to minimize two errors: the brightness error, summed over 
all grid points, and the error in the loop integrals, summed over all elementary square 
paths constructed by connecting the centers of neighboring picture cells 9 . 

The total contribution of the first error term clearly depends on the number of nodes 
in the grid, that is, it depends inversely on e 2 for a fixed image size. We show later 
that the second term, on the other hand, varies directly as e 1 . To make the relative 
contribution of the two terms independent of the grid spacing, we multiply the first term 

7 For a discussion of this issue, within the context of the shape-from-shading problem, see Brooks (1979). 

8 To he exact, it equals zero when z can he written as a polynomial containing only terms of the form 
x'y 3 , for i < 2 and j < 2, for i = j , for i = 0 with j arbitrary, and for j = 0 with i arbitrary. 

9 Strut actually counted each loop integral four times. 



2. Previous work 


7 


by e 2 and divide the second term by e 2 . The quantity to be minimized then becomes 10 

n m n-lm-i 

*=1 3 = 1 t=l j=l 

Here A is a factor that weights the relative contributions of the brightness error and 
the errors in the elementary loop integrals. It can be made small when the irradiance 
measurements are accurate, and the reflectance map is known with precision 11 . 

For the composite error term to be a minimum, the derivatives of the error sum with 
respect to p^i and q^i must be zero. Now, p^i and q y occur in the expressions for e^i, 
e k-l,l, ek-ij- 1 and ekji-i. So performing the indicated differentiations and equating the 
results to zero, one obtains, for .1 < k < n and 1 < l < m, 

e “{Eki - R{pkhgki))Rp(pkhQki) ~ Y e l ek J + e k-i,i - e k-i,i -1 - e k,i-i] — o, 

e (Em — R{pki,qki))R q {pki,qki) — ^ [ e k-i,i + e k-i,i-i — ^k,i -1 — e*,/] = 0, 

where R p and R q are the partial derivatives of R(p, q ) with respect to p and q, as before 12 . 
We can change dummy variables again and gather terms in a particular way to obtain 

( Eij ~~ R[piji Qij)) RpijPij j Qij) T ^\ v ij ~ Rj) ~ 0, 

1 ~ R(Pij> Qij)) R<l(Pij> Qij) + A[/lj-y — ppj] = 0, 

where 

Pij = 4 (Pi+lj+1 ~ PtS&j+l +Pi- lj-1 -Pi+l,j-l) , 

%' ~ i (ft+ij+i “ 9t'-ij'+i A- ft-i,y-i - ft+ij-i) > 
are discrete estimates of the second partial cross derivatives p xy and q xy (times e 2 ), while 

"% “ 4 [(Pt+lj'+l ~ %Pi+lj + Pi+l,j-l) 

T 2(Ptj+i ~ 2 pij + Pi,j- 1) 

T (Pi-lj+l — 2 Pi-ij + Pi—l,y—i)] j 
hij — \ [(ft+lj+l - 2g t , J+ i + qi-ij+i) 

d" 2(ft+ij 2ft y T ft -l.j) 

+ (ft+l,j-l ~ 2ft,y_i + ft_i,y-i)], 

arc discrete estimates of the second partial derivatives p yy and q xx respectively (times e 2 ). 
(Note again that the subscript i in the discrete version corresponds to x in the continuous 
case, while the subscript j corresponds to y.) Strut wrote his result in terms of various 

I'or simplicity, wo assume a rectangular image region liere. There is no loss of generality, however, since 
the sums can be taken over whatever region is desired. 

11 Note, however, that iterative schemes derived directly from the above will become unstable if A is made 
small enough. 

12 Some of the terms are omitted when a point on the boundary is being considered. 
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intermediate expressions, so the equivalence to discrete estimates of partial derivatives 
was not apparent. 

At this point we can isolate the terms in p VJ from one equation, and the term in 
from the other, if we let 

Vij Pij Pi] and h{j — q x j — , 

where p tJ and q x; are given by 

4 [(P»+lj + l — + Pt+lj'-l) + 2(Pi'j + l + Pij- 1) + (Pt-lj + l ~ 2pi-ij + Pi-lj'-l)] , 

\ [(<7t+lj'+l ~ 2<7 jj+1 + Qi-l,j+ 1 ) + Hm+lj + <7+ (ft+ij-i - 2g z j_i + <7;-i,y-i)] > 
respectively. In this way, we obtain 

e 2 

Pij Pij ~ Qij T (^ij R{Piji Qij)) RpiPiji Qij)> 
e 2 

9tj Pij "h i^i] R{Piji qij}') Rq{Pij> Qij)' 

An iterative scheme can now be developed in which the terms p i? and q t j on the left-hand 
side of the equations are considered to be new values that are to be computed by inserting 
the current values into the right-hand sides. Then we obtain: 

| Pij 1 = Pij - tij + J (% - «&&?£■)) R,.(p‘y. A). 

14 +1 == A - +j (% - *(4 4)) *.(44)- 

This scheme appears to work reasonably well, having good stability and convergence 
properties. 

It is clear that one lias to do something special about the boundary, since the above 
result applies only for 1 < k < n and 1 < l < m. On the boundary, different expressions 
apply, which can be obtained by carefully determining which of the terms are missing 
from the result of the initial differentiation. Put another way, the expressions for p;y, cf{ y, 
V{j, and h-ij require the old values of p tJ and q XJ at picture cells bordering on the region 
in which one is applying the iterative scheme. That is, before the scheme can be applied, 
p and q must be known on a border that is one picture cell wide. 

Note that one cannot incorporate occluding boundary information in this scheme 
because, on the occluding boundary, at least one of p and q becomes unbounded. Strat, in 
fact, was forced in his examples to specify the gradient along some closed curve other than 
the occluding boundary. This kind of information is not usually available in applications 
of machine vision. 

2.2. The method of Ikeuchi-Horn 

Ikcuchi and Horn (.1981) were the first to apply the calculus of variations to the shape- 
from-shading problem. They effectively solved a functional minimization problem in 
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recovering object surface orientation. It is known that the occluding boundary provides 
important constraints on the solution of the shape-from-shading problem (I3russ, 1983). 
The difficulty with using the gradient to specify surface orientation is that, as already 
mentioned, at least one of p and q is unbounded on the occluding boundary. 

This problem can be overcome by specifying surface orientation in another way. Con¬ 
sider the mapping from pq space to fg space specified by the equations 

, 2 P , 2( 1 

f — -— 7 ======== === and g = --p====. 

1 + \/T + p 2 + q 2 1 + \/l + p 2 + q 2 

It is easy to verify that / 2 + g 2 < 4 for all visible parts of a surface. The orientation 
of a point on the occluding boundary corresponds to a point on a circle of radius two 
in fg space. Thus occluding boundaries present no difficulties now. The correspondence 
between pq space and the Gaussian sphere of possible orientations can be rationalized 
in terms of the gnomonic projection from the center of the sphere onto a tangent plane. 
Likewise, the correspondence between fg space and the Gaussian sphere can be thought 
of in terms of the stereographic projection from a point on the sphere onto a plane tangent 
to the sphere at the opposite point (see Ikeuchi and Horn, 1981). 

We now seek appropriate / and g values at each point in the image. This we may 
regard as a search for two functions, / and g, defined over fi, that correspond to a smooth 
surface satisfying the image irradiance equation 


E(x,y) = R(f(x,y),g(x,y)). 


(Note that the reflectance map here has been parameterized on / and g.) 

We now develop an appropriate functional. Noting that / and g should ideally cor¬ 
respond to a surface that would produce the image if illuminated the same way as the 
actual surface, we adopt the integral of the brightness error 

E{x,y) - R(f[x,y),g{x,y))^ dxdy. 

We could, at this point, try to add a term that depends on the loop integrals, as Strat 
did. A problem with the use of stereographic coordinates is that the expression for the 
loop integrals becomes complicated. We have 



4/ 

p 4 - f 2 - g 2 


and 


q = 


4.9 

4-/ 2 -9 2 ’ 


so that p y — q x — 0 yields 

fy( 4 + / 2 - g 2 ) - £ 7*(4 - / 2 f g 2 ) + 2 (g y - f x )fg _ 

(4 - / 2 - 

This expression, even when multiplied by (4 — f 2 — g 2 ) 2 , is quite complex and leads to 
even more complicated Euler equations. 

Yet without additional constraint the problem is not well posed. As we saw earlier, 
the minimization of the total brightness error alone does not constitute a well-posed 
problem. In the above case we can choose, at each point (x,y), any / and g for which 
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R(f,g) = E(x,y). In general, there is a one-dimensional family of possibilities—contours 
of constant 11 in fg space. 

We would expect, however, that neighboring points have similar orientations, so that 
a typical “solution” of this form would not be reasonable. Ikeuchi and Horn decided to 
add the measure of “lack of smoothness” given by 

/if ^ + ^ + + 9 ^ ^ dy ‘ 

A solution that produces a small value will be one that keeps the fluctuations in / and g 
small. Adding this term to the brightness error, we obtain the functional 

/if i E ( X,y ^ ~ A (fx + fy + si + S 2 y) dxdy 

that is to be minimized by choosing / and g. Here, again, A is a scalar that assigns a 
relative weighting to the terms. 

The additional expression can be thought of as a regularization term 13 . Such a term 
can be added to a functional in order to obtain a solution in the case that a minimization 
problem does not have a unique solution. 

The Euler equations for this minimization problem can be simplified to read 

(. E-R)R f + XV 2 f = 0 , 

(E - R)Rg + A V 2 .? = 0, 

where Rj and R g are the partial derivatives of R(f,g) with respect to / and g and 

2 d 2 d 2 

dx 2 dy 2 

is the Laplacian operator. 

These Euler equations do not have a unique solution without additional constraint. 
The constraints available to us here are the values of / and g on the occluding boundary. 
Finding the solution of the Euler equations with this particular set of boundary conditions 
usually constitutes a well-posed problem—although this depends on the exact nature of 
the reflectance map, R , and the image, E. 

At this point we introduce a discrete approximation of the Laplacian. The Laplacian 
of a function at a given point is approximately equal to a constant times the difference 
between a local average of the function and its value at the point. The factor of propor¬ 
tionality depends on the way in which the local average is computed. So, for example, if 
we use the simple finite-difference approximation 

~ ^2 [(/*j+1 fi+ij "h fij -1 "b ~ , 

then 

{V 2 /}., * -(/„■ - 


13 Ikeuchi and Horn, however, did not think of the extra term as a regularization term. 
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where the local average, f is given by 

fij ~ l [fij +1 + fi+l,j + fij -1 + • 

The same can be done for g, of course 14 . Using these finite-difference approximations in 
the Euler equations derived above, we obtain 

€ 2 

fij — fij + ^ ~ R{fij,9ij))R/{fij,9ij), 

€ 2 

9ij — 9ij + ^ (Rtj ~~ R[fiji9ij))Rg{fiji9ij)- 

where we have isolated the terms in f tJ and g, j. An iterative scheme can now be developed 
in which these particular terms are considered to be new values to be computed by 
inserting the current values into the remainder of the expression. In this fashion, we 
finally arrive at the scheme 

4 +1 = fa+~ (e,, - r(4 4-)) , 4), 

,2 

. <>T =4+jjf(% -- 

Here, as before, e denotes the spacing between picture cells, while / and g are the local 
averages of / and g. 

This scheme appears to work reasonably well, having good stability and convergence 
properties. We shall see later, however, that the solutions for surface orientation may not 
correspond to an underlying smooth surface and that solutions may be distorted by the 
presence of the regularizing term. The degree of distortion depends on the parameter A. 

2.3. Smith’s approach 

Smith’s method (Smith 1982) was derived by application of the standard calculus to the 
discrete domain. Wc now rationalize his method using the variational calculus. Surface 
orientation can be parameterized in Im space, where 

P 9 

l —- 7 ========== and m = -■=========. 

y/l+p 2 + q 2 y/i + P 2 + <? 2 

This corresponds to an orthographic projection of the Gaussian sphere onto a plane 
tangent to the sphere at one of the poles. We next adopt a regularizing term and minimize 
the functional 

(E(x,y) - R(l(x,y),m(x,y))) 2 + X ((V 2 /) 2 + (V 2 m) 2 ) dx dy. 


14 Slightly bettor results can be obtained using the nine-point approximation of the Laplacian rather than 
the five-point approximation shown here. 
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From the associated Euler equations, wc obtain 

(E - R)Ri - AV 4 / st 0, 

[E - R)R m - A V 4 m - 0, 

where V 4 is the biharmonic operator 15 . 

We need, once again, to impose boundary conditions to avoid ambiguity in the so¬ 
lution. For the biharmonic equation we need to specify l and m on the boundary, as 
well as the normal derivatives of l and m. The normal derivative is the derivative in the 
direction of the outward normal to the boundary curve dfl. Note that while the values of 
l and m on the occluding boundary are known, it may not be obvious what the normal 
derivatives of l and m ought to be. Since they are not specified they must obey the 
appropriate natural boundary condition. 

We can now use the simple finite-difference approximation 

{V 4 l}y « p fry - ty], 

where 

4? — 20 [8(4+1,j T 4-1,j + k,j +1 + hj- 1) 

- 2(4+i, j+i + 4-i.j-i + k-ij+i + 4+ij-i) 

~ (4+2j + 4-2./ + 4,.f+2 + 4,y— 2 )] • 

The same can be done for m, of course 10 . Isolating the terms in 4y and m,+, wc obtain 

e 4 

4j = 4 j + 20 A ~ ^(4j) m ij))Ri{ky< m ij) 1 

e 4 

m ij ~ m ij “h 2QA ~ R(k]i m ij))Rm{kj, m ij)i 

wliich leads to the iterative scheme 

<‘ +1 = ii + ~(E,j - <■))*<(£,<), 

. 1 = < + M 0 % - • <■)) 

The biharmonic equation and its variants are known to require careful treatment. In 
the iterative scheme as written above, for example, the new values are based directly 
on old values. This is called the Jacobi method, and is appropriate for parallel imple¬ 
mentation. But this method is unstable for computational molecules that are discrete 
approximations of the biharmonic operator. A stable iteration can be achieved if one uses 
instead the Gauss-Seidel method, in which the computation of a new value at one picture 

15 The Inharmonic operator can he defined in terms of the Lapladan operator by V 4 (/) = W(V 2 (/)). 

10 Slightly better results can be obtained vising the twentyfive-point approximation of the Inharmonic 
operator instead of the thirteen-point approximation shown here. 
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cell uses the new values of those picture cells already visited in a raster scan of the image. 
This method, however, does not, lend itself to parallel implementation. An alternative 
stabilizing technique depends on the use of smoothing between steps of a Jacobi iteration. 

The more complex boundary conditions mentioned above are reflected in the fact 
that the computational molecule used as the discrete approximation of the biharmonic 
operator requires values for l and m in a band two picture cells wide bordering on the 
region in which the iterative scheme is applied. It is not enough to know the values of l 
and m on the occluding boundary 17 . 

Smith reported difficulties with the above scheme and incorrectly concluded that 
smoothness constraints fail to propagate boundary conditions by more than a few pixels 
in the image. In fact, by a suitable application of the aforementioned stabilization tech¬ 
niques, the scheme can be made to work. Note that fewer problems are encountered if 
[l~ j: + ly T rr? x +rriy) is used in the above functional as the regularization term. This is, in 
part, because the Euler equations then contain the Laplacian operator, for which simple 
iterative schemes exist that are well behaved; but mainly because the treatment of the 
boundary is simpler. 


2.4. Depth from gradient 


A use of the variational calculus in a subsiduary problem arises in the problem of recov¬ 
ering depth from the surface gradient. Let us suppose that we have determined surface 
orientation over the region fi. The relative depth of surface points may be determined 
from the gradient (p, q) by means of the equality 

6 z — p6x -}- q6y, 

that relates infinitesimal changes in x, y and z. Integrating along a curve C from (x 0 ,y 0 ) 
to (x, y), we obtain 

z(x,y) = z(x 0 ,y 0 ) + / (;pdx + qdy). 

JC 

This simple method of integration performs badly when the data are noisy. A depth value 
obtained at some point will, in these circumstances, depend on the integration path that 
was taken to get there. 

It is better to find a best-fit surface z to the given components of the gradient, p and 
q. This we can accomplish by minimizing the functional 


a 


[Zi - p) 2 + (*y - qfdxdy , 


wdiose Euler equation reduces to 


V 2 z = p x + q y . 


Once again, note that this equation does not uniquely specify a solution without further 
constraint. In fact, we can add any harmonic function 18 to a solution to obtain a different 


17 For further details, see for example the discussion of molecular inhibition by Terzopoulos (1983). 

18 A harmonic function satisfies Laplaces equation, V 2 z = 0. 
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solution also satisfying the given Euler equation. In the case here there are no a priori 
boundary conditions given to us. That is, the function sought is not restrained on the 
boundary. The calculus of variations provides us in this situation with natural boundary 
conditions that must be satisfied by the solution. For this particular problem, the natural 
boundary conditions turn out to be (see Appendix) 


(z x ,z y ) ■ n = (p,q) n, 

where 

dy dx 
ds ’ ds 

is a normal vector to the boundary curve <90 and s is arc-length along the boundary. So 
the component of (z x ,z y ) normal to the chosen boundary curve must match the normal 
component of ( p,q ) 19 . 

With these boundary conditions, the solution is still not quite unique, since an ar¬ 
bitrary constant can be added to z without changing the functional. This reflects the 
fact that one cannot recover absolute depth from the gradient (and thus from shading 
information). To get a particular answer, one can fix one of the depth values, or fix their 
average. 

Using the discrete approximation to the Laplacian employed earlier, we obtain the 
iterative scheme 

4 +1 = 4 _ foj + 

where 

z ij — 4 i z i+l,j + z i-lJ T z iJ+l + z i,j- l)j 



is a local average of z, while 


hij — 2 {Pi+ l,y an< ^ v ij ~ 2 QiJ- 1)> 

are estimates of the partial derivatives p x and q y respectively. This is as derived by Horn 
and reported by Ikeuchi (1983). (Note again that the subscript i in the discrete version 
corresponds to x in the continuous case, while the subscript j corresponds to y .) 

In addition to finding the discrete approximation of the Euler equation, we also must 
find the discrete approximation of the boundary condition. This can be done easily, 
provided that the boundary curve is polygonal, with horizontal and vertical segments 
only. This restriction does not provide a problem in our simple situation. Now z x — p 
on vertical segments of the boundary, while z y — q on the horizontal segments. These 
conditions may be translated into 

z i~lJ) — Piji 


19 Note that here s x and z y are the derivatives of the purported solution z(x, y), while p and q are the given 
surface orientation data only with perfect data would these be the same. 
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respectively. These relationships can he used to modify the computation of the average, 
z v j, for points on the edge of the region in which depth is to be reconstructed. Alterna¬ 
tively, these equations can be used to provide phantom depth values on a border of one 
picture cell width around that region. In this case the computation of the average can 
proceed in the same fashion for all points. 


3. Smoothness and integrability 


Methods that attempt to recover shape information encoded in an image usually confine 
their attention to smooth , or piece-wise smooth , solutions. Smoothness, however, is a 
loose term that may be interpreted in many ways. To be specific, we here define a graph, 
z(x,y), to be smooth over a region 0 in the rry-plane if p y = q x , that is, if 


d 2 z _ d 2 z 
dxdy dydx 


V(z,y) <E 0. 


This is a property of C 2 surfaces 20 . Because they must be twicc-differentiable under this 
definition, surfaces that have edges (like polyhedra) are excluded, and it may be argued 
that this accords with our intuitions on smoothness 21 . 

Let us now look more closely at the “lack-of-smoothness” term used in the Ikeuchi- 
Horn method. Suppose that we present a shape-from-shading problem to the program by 
providing it with an image, a reflectance map, and the occluding boundary. The image 
just happens to be that of a Lambertian sphere illuminated by an overhead point source 
at the viewer. This is a well-posed problem with two solutions, a concave bowl and a 
convex ball. It turns out that the algorithm will converge to a somewhat flattened sphere, 
given a planar initial estimate. Interestingly, it converges to almost the same solution 
when given the correct shape initially. That is, the algorithm moves away from the right 
answer. It is interesting to consider why this should be so. 

Recall the functional that is used to derive the Ikcuchi-IIorn method. We are required 
to minimize 

JJ^E(x,y) - R(f( x , y ),g(x,y))'] +A (fx+fy+yl + gl) dxdy. 

It is clear that minimizing the integral of ( E — R(f, y)) 2 is desirable; we wish to make 
the brightness error as small as possible. However, it is not obvious just what is achieved 
by minimizing the remainder of the overall integral. Certainly, it is not smoothness as 


20 For a discussion of this issue, within the context of the shape-from-shading problem, see Woodham 
(1981). 

21 Actually, one may claim that our intuition of smoothness accords better with a less restrictive definition. 
Consider, for example, adjoining a planar surface and a portion of a cylinder along a generator of the 
cylinder. This can be done in such a way that there is no discontinuity in surface orientation, and the 
surface may be considered to be “smooth” (Ikeuchi & Horn, 1981). The second derivatives, however, are 
discontinuous at the join. This suggest that we ought to identify our intuition of smoothness not with 
C 2 , but PC 2 , the set of surfaces that have piece-wise continuous second derivatives. In fact, some may 
even argue that surfaces in 6' 1 , that is, those with continuous first derivatives, are “smooth.” We ignore 
this subtle point in what follows, and restrict attention to surfaces in C 2 . 
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defined above. In fact, if / and g are solutions to the Euler equations for this problem, 
it will in general be the case that there exists no physical surface corresponding exactly 
to the surface orientation specified by / and g (Brooks, 1982). 

The expression (fx+fy+i 7x + 9 y) is instead best regarded as a regularizing term 
that is primarily intended as a means of finding a particularly smooth shape that is close 
to a solution of the original problem (Poggio & Torre, 1984). Different surfaces will give 
rise to different values for 

IL + ^ dx dy - 

Those that fluctuate in depth the least will likely give rise to small values. When a 
shape-from-shading problem is highly ambiguous, in that there is an infinite number of 
possible solutions, a regularizing term is precisely what is needed to get close to one of 
them. If, however, the problem is unambiguous, regularization will usually result in loss 
of accuracy, as the correct solution is unlikely to minimize the integral of the regularizing 
term. 

The distortion due to regularization depends on the parameter A. A large value of A, 
appropriate when the image data is very noisy, leads to large errors, since the emphasis 
will be on producing as smooth a surface as possible, while permitting considerable error 
in brightness. Conversely, a small value for A causes brightness errors to be weighted 
more. In this case, a more undulating surface is acceptable since the contribution of the 
regularizing term to the overall functional is relatively small 22 . 


4. Imposing integrability as a constraint 


In any case, it is desirable to have a shape-from-shading method that neither moves away 
from correct solutions, nor converges to surfaces that are not solutions. To derive such 
a method, we need to impose the smoothness condition defined earlier, instead of using 
regularization. We first consider forcing the solution to satisfy the condition exactly. 

Let us suppose that a shape-from-shading method recovers smooth functions p(x,y) 
and q(x,y ) defined over the image, thereby specifying the gradient. In general, there will 
be no smooth surface that corresponds to this gradient. This is because the functions 
p and q must be related in a special way if they are to correspond to a smooth surface 
(Brooks, 1982). Noting once again that 

d z 3 z 

P{^,y) = g^(^,y) and q(x,y) =—(x,y), 


it follows that, for our earlier definition of smoothness to be satisfied, we must have 


dp 

dy 


{x,y) 


9q 

dx 


(*»!/)» 


or z xy = z yx . This is known as the constraint of integrability. If the gradient does not 
possess this property, there exists no C 2 surface that could give rise to it. Thus we shall 
now attempt to ensure that solutions are integrable. 


22 The iterative scheme becomes unstable, however, when the value of A is reduced too much. 
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4.1. Direct recovery of relative depth 


With the exception of the method of characteristic strips (Horn, 1975), till shape-from- 
shading programs have recovered surface orientation in a separate step, prior to recovering 
relative depth. We stiw earlier an iterative scheme that determines depth values from the 
surface gradient. Of interest here is the direct recovery of depth information, achieved 
without the explicit manipulation of surface orientation. 

Following the guidelines listed earlier, our initial task is to formulate an appropriate 
functional. The brightness error is readily expressed as 


E(x,y) - R(z x (x,rj),z y (x,y)) S j dxdy. 

Now we are to ensure the satisfaction of the constraint z xy = z yx . We might therefore 
consider adding the functional 




dx dy. 


Tt is easy to veryfy that such a term makes no contribution to the subsequent Euler 
equation. This is because the integrand is a divergence expression (Courant Si Hilbert, 
19 5 3) 23 . In our terms, by definition, we seek a smooth surface satisfying the partial 
differential equation. The integrability constraint is redundant. Put yet another way, 
we cannot avoid imposing the integrability constraint if we'look for a scheme that gives 
us z[x i y) directly. This was not the case when we used p and q as parameters. The 
functions p and q had to be related in a special way to satisfy integrability. 

After simplification and reordering of terms, the Euler equation for the brightness- 
error functional alone is 


{RpZ xx T 2 RpR ( jZ xy T R.qZyy) — ( E x Rp T E y R y ) 

— (JF — R){R pp z xx T ‘lRpqZ xy T Rq y z yy ^ } 

where we have used the condition z yx — z xy . Note that p and q replace z x and z y as 
subscripts of R to improve readability. A solution to this equation will give the functional 
an extremal value. By converting the Euler equation to discrete form, employing discrete 
approximations of the derivatives of z, and isolating terms in zij, we obtain the complex 
scheme 


1+1 (Rl + R] -(E- R)(R„ + R„)) 

-k 


ij 


K,(Rl -(E- R)R„) + S| (R„R, - (E - R)R„) + 5 - [E - R)R,„) 

e 2 

T i^xRp + EyRq ) . 


where 


%ij — 4 (■Afl jfl T z i-l.j-l z i+l,j—l)i 


23 A divergence expression in the functional does, however, affect the natural boundary conditions. 
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is a discrete estimate of the cross derivative of z (times e 2 ), and 

K'J — 2^' 2? + 1j ~h l,y)r and Vrj — ^(^ij+i d~ z i,j- i)> 

are horizontal and vertical averages of z respectively. This scheme, unfortunately, is not 
convergent. Other schemes tried also failed. We found little in the literature about how 
one might discover successful iterative schemes for complicated non-linear equations such 
as the one above. Certainly, as far as the variational approach is concerned, the above 
Euler equation must be regarded as fundamental to the problem: the original functional 
is not easily formulated in a more basic way. 

4.2. An alternative approach 

Not surprisingly, if we parameterize the surface on p and q, and impose the integrability 
condition p y = q x , we obtain an Euler equation identical to the one obtained above. The 
functional to be minimized is in this case 

E(x,y ) - R(p(x,y),q(x,y))^) + p(x, y) (p y - q x ) dx dy, 

where p is a Lagrangian multiplier used to enforce the constraint p y = q x . The associated 
Euler equations lead to 

(E — R)R p + ^p y = 0, and ( E — R)R q — \p x = 0. 

In order to eliminate p, we take the (total) derivative of the first equation with respect 
to x and the (total) derivative of the second with respect to y. Adding the results we 
obtain 

{ R pPx + RpRq {Py + Qx) + R qQy) ~ {Ef.Rp + EyRq) 

— ( E R)(R.ppp x T R V q(p y T q x ) T EipiQy) • 

Taken together with the constraint p y = q x , this is the same result as that obtained in 
the previous section. 



5. An integrability penalty term 


It appears to be difficult to extract convergent iterative schemes from Euler equations ob¬ 
tained through the imposition of integrability. Consequently, we now assess the usefulness 
of the penalty term, (p y - q x ) 2 , appearing in the functional 


ff ( E ( x ' y ) " R (P’d)) 2 + A {Py - Qxf dxdy. 


This has the desirable property that if smooth functions p(x,y ) and q(x,y) arc found that 
cause this integral to evaluate to zero, we will, by definition, have solved our problem, 
for the surface will generate the image, and will be smooth everywhere 24 . 


24 Note, however, that we now admit the possibility that p y and q z may be only approximately equal over 
the region 0. 
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The Euler equations for this problem yield 

{E R'jRp 4" A( Pyy — Qxy) — 0 , 

[E - R)R q + A [q xx - p yx ) = 0 . 

Upon isolation of the center term in the discrete approximation of the highest-order, even 
partial derivatives, we arrive at the iterative scheme 



where 

Pii f Pij'—i) and g.-j - o (9i+ l .j T 9i— 1 , 7 ) 

are the vertical average of p and the horizontal average of q, respectively, while p tJ and 
qij are estimates of the cross derivatives (times e 2 ) obtained using the approximations 

Pij — 4 (p*+l,y+l Pi—lJ— 1 ~ Pi— 1 +1 ~~ Pi+lj— i)i 
Qij = +9i-lj-l ~~ 9i —1,7 + 1 — 92 + 1,7-1)5 

respectively. 

This iterative scheme appears to work well. Only very small departures from the 
correct initial solutions have been observed, these being due to the fact that the finite- 
difference expressions are approximations to derivatives. The scheme does not converge 
to a flattened surface as is the case with the Ikeuchi-IIorn method. Rather, we obtain 
asymptotic convergence to the correct solution. Note once again, however, that this 
method requires that the gradient (p,q) be supplied on some closed curve other than the 
occluding boundary. 

This iterative scheme produced very accurate results in tests conducted on synthetic 
images, although, like most shape-from-shading methods, it typically takes many iter¬ 
ations to converge. The observed slow convergence could be alleviated by the recently 
popularized multi-grid technique of processing images and gradient fields at various res¬ 
olutions (Terzopoulos 1984). 

It appears that the use of a penalty term based on a constraint leads to iterative 
schemes that adjust the present estimates in the direction that reduces the penalty term. 
This is in distinction to the behaviour of the schemes that result from attempts to strictly 
enforce the constraint itself. The use of the penalty term gives a scheme some direction¬ 
ality or “push” towards the desired solution. This may be why we were unsuccessful 
in discovering convergent iterative schemes based on the Eider equation derived in the 
previous section. 

5.1. Relationship to Strat’s scheme 

It is interesting to observe how similar the iterative method we derived here is to that 
obtained by Strat. We can sec in retrospect why this should be so, by applying Gauss’s 
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integral formula to Strat’s elementary loop integrals. We have 



(p(x,y) dx + q(x,y)dy) 




dx dy, 


for a simply connected region 72, where the boundary dR is traversed in a counter¬ 
clockwise direction. 

Now, if c is constant in the region 72, then 



where A(R) is the area of the region R. For a smooth surface, p tJ and q x arc continuous, 
so that, for a small enough region R, we can consider them to be nearly constant. That 
is, 


j^(p{x,y)dx + q(x,y)dy)^ = (^J Jjp y - q x ) dx dyj « A{R) Jj (p y - q x ) 2 dx dy . 


So 


^ ,2 
e iJ ~ c 


/ / (Py ( h)~ 

J J6R 


dx dy, 


where <572 is a square region with sides of length e. Consequently, we can consider the 
sum of the error terms squared, 


r, n— 1 m —1 

0 £ £ Ip i.j + Pi+lJ f Qi+l,j + Pi+1,5+1 ~ Pi+1,5+1 — Pt',y+1 ~ Qi,j+l ~ Qi,j]“ > 

1=1 j=1 


or, written more suggestively, 

9 n-l m-1 

(2) X/ [(P»,J + 1 ~ Pij) + (Pi + lJ+l — Pi+l,j ) — (tfi-f l,j ~~ Qi.j) ~ {Qi+lJ+l — QiJ+l)} , 

t-1 J-l 


to be a discrete approximation of 


IIn ^ Py ~ q '^‘ 


dx dy. 


Our final result in the previous section looks a little different from that of Strat, in part 
because we end up with simpler estimates for the second partial derivatives p yy and q xx . 


5.2. Constraints and penalty terms 

We have two equalities: the image irradiance equation, E = 72, and the integrability 
condition, p y = q x . If we enforce both strictly, we obtain Horn’s original characteristic 
strip equations. We have seen that a convergent iterative scheme can be obtained if we 
instead build a functional based on the penalty terms, (E - R) 2 and (p y - q x ) 2 . We also 
described our lack of success in deriving schemes for minimizing the integral of [E — R ) 2 
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while enforcing the constraint p y — q x . We have not yet explored the fourth alternative 
of minimizing the integral of (p y — q x ) 2 while enforcing the constraint E = R. That is, 
minimizing 

JJ^{p y { x , y ) ~ q x { x , y )) 2 + Pi x , y ) (E{x,y) - R(p(x,y),q(x,y))) dxdy. 

The resulting Euler equations are 

{Pyy Pxy] T p{E lV}Rp — 0, 

{Qxx - Pyx ) + p{E - R)R q = 0, 
which, upon elimination of p lead to 

( Pyy Qxi/)Rp — {Rxx Pyx)Rq• 

This equation is to be solved subject to the constraint E — R, of course. We were unable 
to convince ourselves of the utility of pursuit of this particular approach, since we know 
that brightness measurements will be corrupted by noise in practice. 


8. Incorporating occluding boundary information 

One problem not easily coped with is that of dealing with the occluding boundary. Recall 
that the Ikeuchi-Horn method placed considerable emphasis on the ability to be able 
to handle the occluding boundary. So, although we have taken a step forward in the 
above by incorporating integrability, we have also taken a step backwards in that we 
are no longer able to use the occluding boundary. Note, however, that the integrability 
constraint can be expressed using parameterizations that permit incorporation of the 
occluding boundary information. 

Suppose that instead of seeking surface orientation parameterized on p(x,y) and 
q(x,y ), we attempt to recover directly a field of unit normal vectors n(x,y). We need to 
express the integrability constraint in terms of the unit normal and its derivatives. Let 
x, y and z denote unit vectors in the x, y and z directions, respectively. We have that 

n x n y 

p = -- and q = ———, 

n • z n • z 

n • ((riy • z)x — ( n y ■ x)z) __ n • (n^ x(xx z)) 

(n • z) 2 (n • z) 2 ’ 

using the identity (c • a)b — (a • b)c — a x (b x c). Noting that x x z = —y we obtain 25 

Py = --[nnj,y]/(n-z) 2 . 

Similarly, 

q x = +[nn s x]/(n • z) 2 . 


so it follows that 

_ _ (n-x)(nj/-z)-(n-z)(n y -x) 
Py ~ (n-z ) 2 


25 Here [a be] denotes tlie vector triple product a • (b X c). 
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We conclude that the constraint ( p y — q x ) = 0 can be written in the form 

^•([nn I x] + [nn s y])=0. 

As it stands, this form of the constraint will lead to numerical problems in the im¬ 
plementation of an iterative scheme, since (n • z) becomes very small near the occluding 
boundary. It makes sense then to use instead a constraint obtained by multiplying the 
one above by (n • z) 2 , giving 

/' = [n n x x] + [nn^yj. 

One could, of course, tackle this problem using other parametrizations for surface 
orientation, such as / and g. We saw earlier that the integrability constraint expressed in 
terms of / and g is quite complex, and the derivation of the corresponding Euler equations 
somewhat tedious. We felt that the compactness of vector notation provided sufficient 
incentive to tackle the problem the way we did. There is an advantage to using / and 
g, however: one can avoid the redundancy inherent in the use of a vector to represent 
surface orientation, a quantity that has only two degrees of freedom. It is this redundancy 
that leads us to consideration of the pseudo-inverse of a matrix later on. 

6.1. Using a penalty term based on I' 

We are to minimize a functional of the form 

E(x,y) - R.(n(x,y})^ + A/' 2 + fj,( x,y } (n 2 - 1) dxdy. 

Here we iise the Lagrangian multiplier, y,, to enforce the constraint n 2 = 1. The corre¬ 
sponding Euler equation can be simplified tG read 

~{E - R)R n + 2A Vl' n -|- im + \(l' x {n x x) + l' y (n x y)) = 0, 

where 

I'n = n x x x + n y x y 

is the derivative of I' with respect to n, while 

I' x = [n n xx x] + [n x n y y] + [n n yx y], 

I' y = [n liyy y] + [n y n x x] 4- [n n xy x], 

ai - e the derivatives of I' with respect to x and y respectively. We can find the Lagrangian 
multiplier jx by taking the dot product of the Euler equation with n, to give 

fi = (E- R)R n ■ n - 2XI' 2 , 

where we use the fact that I^n = /'. We can now eliminate p by substituting back into 
the Euler equation. The result is 

~{E - R)Ri + 2A/'j' + A (7' (n x x) + l' y (n x y)) = 0, 
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where 


Rn = R n - (R n ■ n)n = n x ( R n x n) 


is the component of R n perpendicniar to n and 


j ' = I' n. 

Note that j'-n = 0, since I ' n -n = I'. In fact, each term in the above equation is orthogonal 
to n. This vector equation thus provides only two constraints on n. The necessary third 
constraint is given by n 2 = 1. 

Now let 


Jx - (ns x n y + n x n yx ) • y and J y = (n y x n x + n x n xy ) • x. 


Then 


~ [ n n xx x] + J x and Iy — [n n yy y] + J y , 
and the Euler equation can be rewritten in the form 

-(E - R)R^ + 2X1'}' + A1 = A (n x x)(n x x) r n xx + (n x y)(n x y) r n 

where 

1 = J x (n x x) + J y (n x y). 

So we can write 


yy 


A (M x n xx + Myn yy ) = A1 + 2A/'j' - (E ~ R)R^, 

where 

Mj = (n x x)(n x x) T and M y - (n x y)(n x y) T 

are the so-called dyadic products of the vectors (n xx) and (n xy) with themselves 20 . We 
now have a non-linear second order partial differential equation for the normal n (x,y). 

We can use the following finite-difference approximations for the derivatives that 
appear: 




and 

n y ~ 2 e ( n bi+i n i,y-i) > 


1 , 

Uxy ~ 4? v^+bi+i + n i 

-lj-l 

— n *~i,y+i — , 

as well as 




ftxx 

~ ~2 ( n f+l,J ~ 2n ij + n t-l,y) 

and 

1 , 

n yy ~ ^2 ( n *j+i ~ 2n ij + n i,j 

or 

^ (h„ - nij) 

and 

„ 2 /_ \ 
n yy ~ £ 2 V v *.7 n ij) > 

where 


and 

Vij = 2 ( n ij+l *b n td-t) 


20 These dyadic products are matrices of rank one. 
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are horizontal and vertical averages of n respectively. 

We now develop an iterative scheme based on the isolation of the center term in the 
discrete approximations of the highest-order, even partial derivatives. For convenience, 
let mij, say, be the new value of the normal to be calculated in the iterative step. Then 


2 - 2 

— (h ij - m ij) and n yy « -j (v, 


m tjj 


If we let M = M x + M y , then the new value is obtained using the equation 

2 2 

Mm = (M x h + M y v) - U - e 2 /'j' + ~{E - R)R^ 


and the constraint m 2 = 1. Here we omit subscripts in order to simplify the notation. 
Let r denote the right hand side of the equation above. All terms in r can be easily 
computed using the old estimate of the normal, n, in the expressions for M x , My, h, v, 
1, j', R and R 4 . The remaining problem is the solution of the equation for the new 
estimate of the normal, m. 


6.2. Solving the equations Mm = r and m • m = 1 

The equation Mm = r is underdetermined, since M here only has rank two. The matrix 
is singular and so does not have an inverse in the usual sense. There are an infinite 
number of solutions that can be written in terms of the pseudo-inverse , M + . They are 

m = M + r + (l — M + M) x, 

for arbitrary x (Albert, 1982), where I is the 3x3 identity matrix. Of these solutions, 
we seek the one with unit norm, m 2 = 1. 

The pseudo-inverse of a matrix M can be defined as the limit 

M 4 = lim (m 2 ’M -|- <5 2 l) 1 M t . 

6—>0 V J 

Alternatively, it can be defined using the conditions of Penrose (Albert, 1982), which 
state that the matrix M + is the pseudo-inverse of the matrix M if, and only if, 

® MM + and M + M are symmetric, and 

• M 4 MM + = M, as well as, 

• MM+M = M+. 

The pseudo-inverse may also be found using spectral decomposition. The eigenvectors of 
the pseudo-inverse are the same as those of the original matrix, while the corresponding 
non-zero eigenvalues are the inverses of the non-zero eigenvalues of the original matrix. 
Now, in our case, 

M = (nx x)(n x x) T + (n x y)(n x y) r , 

so one can show that 

1 

Kl ) 2 


M + = I - nn r + 


(n x z)(n x z) T . 



6. Incorporating occluding boundary information 


25 


It is also possible to verify that 


M+M = I - nn r , 


from which it follows that 

I - M+M = nn r , 


and so 


m = M+r -f un , 


for some u chosen to make m 2 = 1 . In our case r ± n, so we can further simplify matters 
by noting that 

M+r = r + -fr n zl (n x z). 

(n • zy 

Let this be called p. Since p i. n, we have that m 2 = p 2 + iA This completes the 
calculation of the new estimate of the normal, m. The only potential problem occurs 
when |M+r| > 1 , as may happen when the current estimate of the solution is far from 
the correct one. In this case it is advisable to limit the adjustment of the local normal 
away from its previous value, n 27 . 


6.3. Using a penalty term based on I 

Implementations of the above iterative scheme work well except for minor problems near 
the occluding boundary. What happens is that the components of n s and n ?/ become 
unbounded on the occluding boundary, so that the individual terms in 

I' - [n n x x] T [n tiy y] . 

tend to become very large 28 . It may be better to use the slightly more complicated 
expression 

/ = (n • z)/' - (n • z) ([nn x x] + [nn 9 y]) . 

This can be viewed as the difference of two quantities that remain bounded, provided 
that the curvature of the surface is bounded. 

W'e now are to minimize a functional of the form 

E(x,y) - R(n(x,y))^j + A/ 2 + p(x,y) (n 2 - 1) dxdy. 

The corresponding Euler equation can be simplified to read 
-{E - R)R n + A/(/'z + 2(n • z )I' n ) + /m + 2A/k 

+ A ( n • 2 ) 2 (4( n x x ) + 4(n x y)) = 0, 

27 Also, note that there are theoretically two solutions for v, one positive ami one negative. The positive 
value leads to a now estimate close to the previous one, while the negative value gives rise to one almost 
opposite to the old one. It is clear that one should use the positive root. 

28 The problem is different, of course, from the one we encountered earlier when using the gradient to 
parameterize surface orientation. The components of the gradient, p and q become infinite on the 
boundary, while n remains finite. 
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where 

k = (n z • z)(n x x) + (n y • z)(n x y), 
and I' n , I' x and /' are as defined before. 

We can find the Lagrangian multiplier /i by taking the dot product of the Euler 
equation with n. Thus we have 

!i = {E~ R)R n -n-3XI 2 . 


Now we eliminate /j by substituting back into the Euler equation. The result is 
~{E - R)Rn + A/j + 2A/k + A(n • z) 2 (4(n x x) + I ' y {n x y)) = 0, 

where 

j = I'i + 2(n • z)4 - 3 In, 

and /?„ is the component of R n perpendicular to n as before. Now I' n n —- I', so j n = 0. 
In fact, each term in the above equation is orthogonal to n. This vector equation thus 
provides only two constraints on n. The necessary third constraint is again given by 
n 2 = 1. 

Now let J x and J y be defined as before. Then the Euler equation can be rewritten in 
the form 


■{E - R)R£ + A/j + 2A/k + A(n • z) 2 l 

= A(n • z) 2 (n x x)(n x x) T n XI + (n x y)(n x y) T n 


yy 


where 1 — J x (n x x) + J y [n x y). So we obtain 


A(n • z) 2 (M x n xx + M y n yy ) ~= A(n ■ z) 2 I + 2A/k + A/j - (E - R)R± t 


where M x and are defined as before. We now have a non-linear second order partial 
differential equation for the normal n(x,y). 

Note that both sides of this equation are orthogonal to n, since 1 • n = 0, k • n - 0, 
j • n = 0, and R^ • n = 0. So the equation provides two constraints only, with the third 
coming from n 2 = 1. 

If we use the same discrete approximations as before, and isolate the central value 
in the finite-difference approximations of the highest order even partial derivatives, we 
obtain 


Mm = (M*fi + M„v) 


-/' k 


,2 


2^ (n-z)*** 2(n • z) J ' 2A(n-z) 2 


I'j + 


:(E - R)Ri. 


We once again obtain an underdetermined equation, of the form M m = r, together with 
a constraint m 2 — 1. We can solve for the new estimate of the surface normal using the 
pseudo-inverse of the matrix M, as before. 

It is curious that several of the terms involve division by (n • z), a term that becomes 
large near the occluding boundary. Wc multiplied the penalty term by this expression 
in the first place in order to avoid problems near the occluding boundary. Apparently, 
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however, the terms so affected are all small near the occluding boundary anyway. In 
fact, we determined experimentally that several of the terms on the right hand side are 
very small compared to the others, particularly as one approaches the correct solution. 
We found that one can leave them out without noticably affecting convergence, or the 
surface arrived at ultimately. Preliminary testing of the scheme on synthetic images 
yielded promising results. Comprehensive assessment of the performance of the two 
schemes has, however, been left for future work. 

7. Summary 

The shape-from-shading problem was regarded here as one of finding a surface that 
minimizes an integral expression involving the brightness error. The expression we used 
has the form of a functional measuring the departure of a hypothesized surface from 
a solution surface. Iterative schemes for solving the shape-from-shading problem were 
based on the appropriate Euler equation. 

We reviewed the use of a regularization term in an existing iterative scheme. Regu¬ 
larization techniques allow one to obtain results when faced with ill-posed problems. We 
argued, however, that the addition of a regularization term is not appropriate when one 
is dealing with a well-posed problem. The additional term tends to flatten and distort 
the solution. 

We next discussed the fact that surface orientation must satisfy an integrability con¬ 
straint if it is to correspond to an underlying smooth surface. The method using the 
regularization term does uot guarantee this. We attempted to use the integrability con¬ 
straint instead of a regularization term, but failed to find convergent iterative schemes 
for solving the resulting Euler equations. 

A convergent iterative scheme was obtained, however, when, instead of enforcing 
integrability, we introduced a penalty term derived from the integrability constraint. It 
seems that the penalty term provides the iterative process with a “sense of direction” that 
helps it head towards the solution. This approach allows one to recover surface gradients 
that are approximately integrable. The scheme so derived was shown to be similar to 
that obtained in the discrete domain by Strat. A drawback of his scheme is its inability 
to incorporate occluding boundary information. 

We overcame this difficulty by employing a different parametrization for surface ori¬ 
entation. The integrability penalty term can bo expressed in terms of the unit surface 
normal and its derivatives. Subsequent application of the variational calculus proved to 
be somewhat involved, but two usable iterative schemes were finally obtained. Initial 
tests indicate that they perform well. Our new schemes are the first to make use of the 
integrability constraint while allowing incorporation of the occluding boundary normals. 
Future work will assess the relative performance of the new schemes in detail. 
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Appendix 


A. Minimizing Functionals 

Let F(x,y,z,z x ,z y ) measure the distance of a surface, 2 , from a satisfactory solution at 
a point ( x,y ). For now, assume that F is dependent not only on 2 , but also on the first 
partial derivatives z x and z y . Given that we seek a surface defined over some region Q in 
the plane, and that F is everywhere non-negative, we may regard 

h{z) = F(x,y,z,z x ,z y ) dxdy 

as an overall measure of error whose value is to be minimized. This is not a conventional 
minimization problem since we search over a space of functions and not a region of 
coordinate space. The value of p depends on the choice of the function z, and for this 
reason p is termed a functional. Minimizing p is a problem in the calculus of variations. 

A fundamental result of the calculus of variations is that the extrema of functionals 
must satisfy an associated Eider equation over the domain of interest. For the above form 
of the functional, the equation is 


F --F 8 F 

z dx Tz * dy lz • 


0 . 


This is a necessary condition for the existence of an extremum, 2 (p. 185, Courant & 
Hilbert, 1953). It is not a sufficient condition. Note that local minima, global minima, 
local maxima, global maxima, and inflexion points are all examples of extrema. 

It will prove useful to note two other Euler equations corresponding to other forms 
for F. In the event that F is dependent also on the second partial derivatives as in 


h{z) = JJ^F(x,y,z,z x , 


*y,z. 


XX ) Z x y, Zyy 


) dx dy, 


the Euler equation expands to 
d 


F-^-F ~—F A-^-F , 82 

Z dx lzx dy Izy + d x 2 Zzx + dxdy 


F z 


d 2 r, 

1 Ay 2 Fz ' 


0. 


Sometimes we seek a surface that is parameterized not in terms of relative depth, but 
in terms of surface normals. Two parameters are needed in this case. If the functions p 
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and q are used to describe surface orientation and if the associated functional incorporates 
their first partial derivatives in x and y, the expression to be minimized then takes the 
form 

h{p,q) = J J^F{x,y,P,(},Px,Py,qx,qy)dxdy, 

which has two corresponding Euler equations given by 


F v - 

d F p 

d 

P 

dx Px 

dy 

F a - 

d 

- F a 

d 


dx qi 

dy 


Py 


1y 


0 , 

0 . 


In general, these constitute a pair of coupled partial differential equations in p and q. A 
pair of functions satisfying these equations will generate an extremum of 1 3 . As before, 
the extremum may be a minimum, maximum or inflexion point. 

Note that if the functional involves higher derivatives of p and q, the Euler equations 
generalize in a straightforward way. Courant and Hilbert (1953) provide an excellent 
chapter on the variational calculus. 


B. Boundary Conditions 

In general, a problem involving partial differential equations is ill posed in the absence 
of suitable boundary conditions because the solution is not unique without additional 
constraint. The type of boundary condition that ensures a given problem is well posed 
depends on the particular type of partial differential equation. There may be more than 
one way of adding boundary conditions to a partial differential equation in order to force 
a unique solution. 

In our case, boundary conditions may be given as part of the basic minimization 
problem. That, is, the solution sought must minimize the functional subject to additional 
constraints, such as prescribed values on the boundary of the region of integration. In 
the case that the function is not constrained on the boundary, however, the calculus of 
variations provides so called natural boundary conditions that the solution must satisfy. 
For example, for the functional I\ given above, a suitable boundary condition is the value 
of z along the boundary of the region H. The natural boundary condition in this case 
is just 

{F Zx ,F Zy )-n = 0, 

where the normal to the boundary, n, is given by 


n = 


dy dx\ 
ds ’ ds ) ’ 


and s is arc-length measured along the boundary <9fi. So the component of the vector 
(F Zx ,F Zy ) normal to the boundary should be everywhere zero. 

In the case of the functional I 3 given above, the values of p and q along the boundary 
will usually be suitable. The natural boundary conditions in this case happen to be 

( F P* > F Py ) • 11 = 0 an(1 { F 'h > F q J • 11 = 0 - 
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C. Regularizing Terms 

At times, the problem of minimizing a functional is not well posed as there is an infinite 
number of solutions, even with constraints on the boundary. One can then find a surface 
that is close to a solution, while minimizing some measure of departure from smoothness, 
by regularization (sec Poggio & Torre, 1984). The regularization method of interest to us 
here involves the addition of a regularizing term to the functional. If we deal with the 
problem of recovering a surface z(x, y ) from shading, we may wish to include a regularizing 
term such as the square Laplacian appearing in 

JJ (VV ) 2 dxdy, 

or the quadratic variation in 



Each of these has the desirable property of rotational invariance (Brady and Horn, 1983). 
The lower order rotationaily-symmetric regularizing term 


IL 


(zl +Zy)dx dy, 


leads to excessive flattening of solutions. The latter form may well be appropriate when 
applied to surface orientation parameters, such as p and q , or / and g, but this is not the 
case with a depth parameter, such as z. 


D. Enforcing Constraints 

Sometimes we seek a minimum of a functional subject to some independent constraint. 
Suppose, for example, that we arc required to minimize the previously defined h(z), 
subject to the constraint that 

g(x,y,z,z x ,z y ) = 0. 

In this case we may use the Lagrangian multiplier method in which we minimize not I\ 
but the augmented functional 

h{z) = JJ^F(x,y,z,z x ,z y ) + p{x,y)g(x,y,z,z x ,z y )dxdy. 

We now seek solutions to the associated Euler equation. Note that the Lagrangian mul¬ 
tiplier p is a function of x and y and must be treated as such when deriving the Euler 
equation. 

If we differentiate the functional with respect to p(x,y), for a particular x and y, and 
set the result equal to zero, we get back the original constraint equation. This equation 
is required to help solve for p, something we typically have to do in order to eliminate it 
from the Euler equation. At times, this may take some skill. More importantly, however, 
the equations that result often do not suggest convergent iterative schemes. In any case, 
a solution of the resulting Euler equation will have the property that g(x,y,z,z x , z y ) = 0, 
with I 4 having an extremal value on the manifold g(x,y,z,z x ,z y ) — 0. 
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E. Penalty terms based on constraints 

In view of the difficulties experienced when attempting to impose constraints exactly, we 
often consider an alternative method. In this approach, a penalty term derived from the 
contraint is employed. Thus we might rely on the Euler equation corresponding to the 
functional 

h{z) = JJ" F(x,y, z, z x , zy) + A [g{x, y, z, z x , z y )] 2 dx dy. 

Here, A is a scalar that aligns the arbitrary scales of F and g. Alternatively, it may be 
regarded as a weighting of the relative importance of the components of the functional. 
It is, of course, not necessary to square g if it is already guaranteed to be non-negative 
over Q for all functions z. 

Solutions to the Euler equation for h 0 now specify surfaces that generate an ex¬ 
tremal value of /q. However, these surfaces will not, in general, satisfy the constraint 
g(x,y,z,z x ,z y ) = 0 exactly. Rather, it will be the case that the value of g is small, 
along with the values of the other expressions being minimized. This is usually an ac¬ 
ceptable compromise. More often than not, this approach proves more tractable than the 
Lagrangian method as there is no multiplier to be eliminated. 



